Analysis_01_prone_incidence_risk_reliability
Chad H Hochberg and Anna Barker
2025-08-22
#Specify Project Location and File Type
project_location <- '~/Library/CloudStorage/OneDrive-JohnsHopkins/Research Projects/CLIF/CLIF_Projects/ProneProjects/ProningEpi/ProneSevereARF_Output'
file_type <- 'csv'
#Create Sub Folders within Project Folder
# Check if the output directory exists; if not, create it
if (!dir.exists(paste0(project_location, "/summary_analysis"))) {
dir.create(paste0(project_location, "/summary_analysis"))
}
set.seed(3221984)df_agg <- read.csv(paste0(project_location, '/summary_output/aggregate_data.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
clif_hospital_description <- read.csv(paste0(project_location, '/clif_hospital_description.csv'))
#Merge Hospital Data Into DF AGG
df_agg <- df_agg |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Do this Before Filtering
df_unique_hosp <- df_agg |> distinct(hospital_id, Hospital_Type)
cat('\nHow Many Academic v Community Hospitals are In THese Data')##
## How Many Academic v Community Hospitals are In THese Datatable(df_unique_hosp$Hospital_Type)##
## Academic Community
## 12 25df_hosptype_summary <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.df_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 13 10 12
## 2 Post-COVID 12 11 11
## 3 Pre-COVID 8 8 10#Filter Out Hospitals with Less Than 10 Patients
df_agg <- df_agg |>
#Filter OUT HOspitals-Periods with < 10 observations
filter(n_patients>=10)
#Now After Filtering
cat('\nAfer Filtering Out Hospitals with >= 10 Patients \nHow Many Academic v Community Hospitals are In THese Data')##
## Afer Filtering Out Hospitals with >= 10 Patients
## How Many Academic v Community Hospitals are In THese Datatable(df_agg$Hospital_Type)##
## Academic Community
## 33 50df_hosptype_summary_post <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.cat('\nProning Counts in Hospitals with >= 10 Patients')##
## Proning Counts in Hospitals with >= 10 Patientsdf_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 12 10 12
## 2 Post-COVID 7 11 11
## 3 Pre-COVID 3 7 10#Academic Size Summary
df_hospsize <- df_agg |>
group_by(Hospital_Type) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description by Hospital Type')##
## Hospital Size Description by Hospital Typedf_hospsize## # A tibble: 2 × 4
## Hospital_Type small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 Academic 0 3 30
## 2 Community 22 25 3#Proning By Hospital Size
df_hospsize <- df_agg |>
group_by(icu_bed_cat, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, icu_bed_cat)## `summarise()` has grouped output by 'icu_bed_cat'. You can override using the
## `.groups` argument.cat('\nProning Counts by Hospital Size and Period')##
## Proning Counts by Hospital Size and Perioddf_hospsize## # A tibble: 9 × 6
## # Groups: icu_bed_cat [3]
## icu_bed_cat study_period n_Hospitals n_patients n_proned percent_proned
## <fct> <chr> <int> <int> <int> <dbl>
## 1 small COVID 12 433 300 69.3
## 2 medium COVID 10 623 283 45.4
## 3 large COVID 12 1906 738 38.7
## 4 small Post-COVID 7 128 53 41.4
## 5 medium Post-COVID 11 335 77 23.0
## 6 large Post-COVID 11 1109 179 16.1
## 7 small Pre-COVID 3 43 6 14.0
## 8 medium Pre-COVID 7 156 23 14.7
## 9 large Pre-COVID 10 979 67 6.84#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned <- df_agg
df_agg_proned <- df_agg_unproned %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned <- df_agg_unproned |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg <- bind_rows(df_agg_proned, df_agg_unproned)
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
global_agg_byoutcome_period <- read.csv(paste0(project_location, '/summary_output/global_aggregate_byoutcome_period.csv')) |>
mutate(prone_log_odds = log(prone_rate_adjust / (1 - prone_rate_adjust))) |>
select(hospital_id, prone_12hour_outcome, prone_log_odds, study_period) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Back with DF Agg
df_agg <- left_join(df_agg, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, study_period) |>
arrange(hospital_id, study_period, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, study_period, prone_12hour_outcome)`#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Also Includes Hospital Type and Hospital Type and Period Interaction
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_full <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_full)
summary(agg_brm_period_full)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.90 0.12 0.69
## sd(study_periodPostMCOVID) 1.01 0.18 0.71
## sd(study_periodPreMCOVID) 0.99 0.26 0.58
## cor(study_periodCOVID,study_periodPostMCOVID) 0.77 0.11 0.51
## cor(study_periodCOVID,study_periodPreMCOVID) 0.60 0.20 0.13
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.52 0.23 0.01
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.17 1.00 1302
## sd(study_periodPostMCOVID) 1.40 1.00 1709
## sd(study_periodPreMCOVID) 1.60 1.00 1942
## cor(study_periodCOVID,study_periodPostMCOVID) 0.93 1.00 1957
## cor(study_periodCOVID,study_periodPreMCOVID) 0.91 1.00 2005
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.89 1.00 1788
## Tail_ESS
## sd(study_periodCOVID) 2288
## sd(study_periodPostMCOVID) 2646
## sd(study_periodPreMCOVID) 2578
## cor(study_periodCOVID,study_periodPostMCOVID) 2703
## cor(study_periodCOVID,study_periodPreMCOVID) 1953
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 2730
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## study_periodCOVID -0.11 0.23 -0.57
## study_periodPostMCOVID -1.20 0.29 -1.79
## study_periodPreMCOVID -2.10 0.33 -2.78
## Hospital_TypeCommunity 0.83 0.30 0.27
## study_periodPostMCOVID:Hospital_TypeCommunity 0.12 0.32 -0.51
## study_periodPreMCOVID:Hospital_TypeCommunity -0.16 0.48 -1.16
## u-95% CI Rhat Bulk_ESS Tail_ESS
## study_periodCOVID 0.34 1.00 596 1192
## study_periodPostMCOVID -0.65 1.00 899 1821
## study_periodPreMCOVID -1.49 1.00 1748 2075
## Hospital_TypeCommunity 1.42 1.00 830 1272
## study_periodPostMCOVID:Hospital_TypeCommunity 0.76 1.00 2021 2992
## study_periodPreMCOVID:Hospital_TypeCommunity 0.78 1.00 2285 2279
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## pling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.689 seconds (Warm-up)
## Chain 2: 3.515 seconds (Sampling)
## Chain 2: 6.204 seconds (Total)
## Chain 2:
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.702 seconds (Warm-up)
## Chain 4: 3.636 seconds (Sampling)
## Chain 4: 6.338 seconds (Total)
## Chain 4:
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.875 seconds (Warm-up)
## Chain 1: 3.676 seconds (Sampling)
## Chain 1: 6.551 seconds (Total)
## Chain 1:#Now Generate Odds Ratios for Communtiy vs Academic Proning by Period
full_posterior <- as_draws_df(agg_brm_period_full) |>
mutate(pre_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPreMCOVID:Hospital_TypeCommunity`),
covid = exp(b_Hospital_TypeCommunity),
post_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPostMCOVID:Hospital_TypeCommunity`))
summary_posterior <- data.frame(
variable = c(
'Pre-COVID Community Hospital',
'COVID Community Hospital',
'Post-COVID Community Hospital'
),
median_or = c(median(full_posterior$pre_covid),
median(full_posterior$covid),
median(full_posterior$post_covid)),
conf.low = c(quantile(full_posterior$pre_covid, p=0.025),
quantile(full_posterior$covid, p=0.025),
quantile(full_posterior$post_covid, p=0.025)),
conf.high = c(quantile(full_posterior$pre_covid, p=0.975),
quantile(full_posterior$covid, p=0.975),
quantile(full_posterior$post_covid, p=0.975)))
print(summary_posterior)
write.csv(summary_posterior, paste0(project_location, '/summary_output/summary_posterior.csv'))####Predictions Command from Marginal Effects with 4000 posterior draws
reliability_adjust_full <- predictions(agg_brm_period_full,
type = 'response', ndraws=4000, re_formula = ~ (0 + study_period | hospital_id))
reliability_adjust_full <- as_tibble(reliability_adjust_full) |>
left_join(df_agg %>%
select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds', 'prone_log_odds', 'study_period', 'prone_12hour_outcome'),
by = c('hospital_id', 'study_period', 'prone_log_odds')) |>
#Need to Recombine Prone_12hour_outcome Estimates into 1 Per Hospital and Period
group_by(hospital_id, study_period) |>
mutate(
estimate_add=sum(estimate, na.rm=T),
ci_low_add=sum(conf.low, na.rm=T),
ci_hi_add=sum(conf.high, na.rm=T),
n_patients=sum(n_patients, na.rm=T)
) |>
arrange(hospital_id, study_period, -prone_12hour_outcome) |> #THis Arranges so Row Number 1 is the prone_12hour_outcome==1
filter(row_number()==1) |>
ungroup() |>
mutate(adjust_rate=estimate_add/n_patients) |>
mutate(ci_low=ci_low_add/n_patients) |>
mutate(ci_hi=ci_hi_add/n_patients) |>
mutate(period_rank=fcase(
study_period=='Pre-COVID', 0,
study_period=='COVID', 1,
study_period=='Post-COVID', 2
)) |>
#Orders for Graph
arrange(period_rank, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
mutate(actual_rate=observed_prone/n_patients) |>
mutate(grouping='Actual Rate') |>
left_join(df_agg %>% select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds')) |>
distinct()## Warning in left_join(as_tibble(reliability_adjust_full), df_agg %>% select("hospital_id", : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 13 of `x` matches multiple rows in `y`.
## ℹ Row 13 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Joining with `by = join_by(hospital_id, Number_of_ICU_Beds)`## Warning in left_join(mutate(mutate(mutate(arrange(mutate(mutate(mutate(mutate(ungroup(filter(arrange(mutate(group_by(left_join(as_tibble(reliability_adjust_full), : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 151 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.break_points <- reliability_adjust_full |>
mutate(pre_covid_num=fcase(
row_number()!=1 & study_period=='Pre-COVID' & lead(study_period)=='COVID', hospital_rank),
post_covid_num=fcase(
study_period=='COVID' & lead(study_period)=='Post-COVID', hospital_rank
)) |>
summarise(
pre_covid_num=mean(pre_covid_num,na.rm=T)+0.5,
post_covid_num=mean(post_covid_num, na.rm=T)+0.5)
caterpillar_hospital_type <- ggplot(reliability_adjust_full,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.1),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals",
y = "Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
geom_vline(xintercept = break_points$pre_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 7, y = 0.95, label = 'Pre Pandemic', fontface = 'bold') +
geom_vline(xintercept = break_points$post_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 35.5, y = 0.95, label = 'Pandemic', fontface = 'bold') +
annotate("text", x = 68, y = 0.95, label = 'Post Pandemic', fontface = 'bold') +
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.5))
caterpillar_hospital_type
ggsave('caterpillar_hospital_byperiod_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
#With Actual Rate
#Add the Actual Value
with_rate <- caterpillar_hospital_type +
geom_point(aes(y=observed_prone/n_patients), color = 'maroon', alpha = 0.25)
with_rate
ggsave('caterpillar_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)#Stack Caterpillar Plots Vertically
vertical_data <- reliability_adjust_full |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
ungroup()
vertical <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.1, 0.90),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 9)) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1,
scales='free_x')## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ggsave('vertical_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
verticald <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.75, 0.85),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 12),
legend.key.height = unit(1.5, "lines"),
legend.key.width = unit(2.2, 'lines')) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1)
ggsave('vertical_caterpillar_option_D.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
stacked_data <- vertical_data |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=if_else(
study_period=='COVID', row_number(), NA_real_
)) |>
group_by(hospital_id) |>
fill(hospital_rank, .direction = 'updown')
stacked <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, color = study_period)) +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= study_period), alpha = 0.4) +
geom_point(aes(color=study_period, shape=Hospital_Type), size = 2.5) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank())
ggsave('stacked_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked2 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, linetype=study_period), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) +
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Demuth") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar2.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked_data <- stacked_data |>
group_by(study_period) |>
mutate(period_n=sum(n_patients)) |>
ungroup()
stacked3 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, alpha=period_n), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) + scale_alpha(guide="none")+
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2.25) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Egypt") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Pandemic Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals",
shape="Study Period",
color= "Hospital Type"
) +
scale_shape_manual(values= c("circle", "triangle", "square"), labels=c("Pandemic", "Post-pandemic", "Pre-pandemic"))+
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar3.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).posterior <- as_draws_df(agg_brm_period_full)
# Extract the standard deviations (SDs) for each study_period random effect
sd_pre <- posterior$`sd_hospital_id__study_periodPreMCOVID`
sd_covid <- posterior$`sd_hospital_id__study_periodCOVID`
sd_post <- posterior$`sd_hospital_id__study_periodPostMCOVID`
#3 Calculate MOR for Each Period
# MOR calculation function
calc_mor <- function(sd_samples) {
exp(sqrt(2) * sd_samples * qnorm(0.75)) # qnorm(0.75) ~ 0.674
}
# Compute MORs
mor_pre <- calc_mor(sd_pre)
mor_covid <- calc_mor(sd_covid)
mor_post <- calc_mor(sd_post)
# Differences Using The Latter Periods as Reference
diff_covid_pre <- mor_covid - mor_pre
diff_post_pre <- mor_post - mor_pre
diff_post_covid <- mor_post - mor_covid
# Posterior probabilities that the Difference in mOR is > 0 (Bayesian 'p-values')
prop_covid_gt_pre <- mean(diff_covid_pre > 0)
prop_post_gt_pre <- mean(diff_post_pre > 0)
prop_post_gt_covid <- mean(diff_post_covid > 0)
# Credible intervals
cr_covid_pre <- quantile(diff_covid_pre, c(0.025, 0.5, 0.975))
cr_post_pre <- quantile(diff_post_pre, c(0.025, 0.5, 0.975))
cr_post_covid <- quantile(diff_post_covid, c(0.025, 0.5, 0.975))
# Summarize
mor_summary <- data.frame(
Period = c("Pre-MCOVID", "COVID", "Post-MCOVID"),
MOR_Mean = c(mean(mor_pre), mean(mor_covid), mean(mor_post)),
MOR_Median = c(median(mor_pre), median(mor_covid), median(mor_post)),
MOR_2.5 = c(quantile(mor_pre, 0.025), quantile(mor_covid, 0.025), quantile(mor_post, 0.025)),
MOR_97.5 = c(quantile(mor_pre, 0.975), quantile(mor_covid, 0.975), quantile(mor_post, 0.975))
)
print(mor_summary)## Period MOR_Mean MOR_Median MOR_2.5 MOR_97.5
## 1 Pre-MCOVID 2.670428 2.497792 1.732122 4.604875
## 2 COVID 2.375937 2.336800 1.927092 3.052762
## 3 Post-MCOVID 2.649824 2.578343 1.968028 3.786942#Summarize the Differences in Median Odds Ratios by Period
diff_mor <- data.frame(
contrast = c('COVID-Pre', 'Post-COVID', 'Post-Pre'),
MOR_diff_mean = c(mean(diff_covid_pre), mean(diff_post_covid), mean(diff_post_pre)),
MOR_diff_median = c(median(diff_covid_pre), median(diff_post_covid), median(diff_post_pre)),
MOR_diff_2.5 = c(quantile(diff_covid_pre, 0.025), quantile(diff_post_covid, 0.025), quantile(diff_post_pre, 0.025)),
MOR_diff_97.5 = c(quantile(diff_covid_pre, 0.975), quantile(diff_post_covid, 0.975), quantile(diff_post_pre, 0.975)),
prob_diff_gt_0 = c(prop_covid_gt_pre, prop_post_gt_covid, prop_post_gt_pre)
)
print(diff_mor)## contrast MOR_diff_mean MOR_diff_median MOR_diff_2.5 MOR_diff_97.5
## 1 COVID-Pre -0.29449118 -0.13969104 -2.2474417 0.7870697
## 2 Post-COVID 0.27388686 0.23031086 -0.5037268 1.3189232
## 3 Post-Pre -0.02060432 0.08311376 -2.0002086 1.4344507
## prob_diff_gt_0
## 1 0.39325
## 2 0.72325
## 3 0.54450#Visualize
mor_df <- data.frame(
MOR = c(mor_covid, mor_pre, mor_post),
Period = factor(rep(c("COVID", "Pre-MCOVID", "Post-MCOVID"), each = length(mor_covid)))
)
ggplot(mor_df, aes(x = MOR, fill = Period)) +
geom_density(alpha = 0.5) +
labs(
title = "Period-Specific Median Odds Ratios (MOR \n (Fully Adjusted))",
x = "Median Odds Ratio (MOR)",
y = "Density"
) +
scale_x_continuous(seq(1,20, by=1), limits =c(1,20),
name='Distribution of Median Odds Ratos by Study Period') +
theme_minimal()
ggsave('mor_distribution_fully.pdf',
path = paste0(project_location, '/summary_output/graphs/'))## Saving 7 x 5 in image#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Adjusted for COVID Status
#First Prepare Dataset
#Open Aggregate Data with COVID
df_agg_covid <- read.csv(paste0(project_location, '/summary_output/aggregate_data_wcovid.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Hospital Data Into DF AGG
df_agg_covid <- df_agg_covid |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned_covid <- df_agg_covid
df_agg_proned_covid <- df_agg_unproned_covid %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned_covid <- df_agg_unproned_covid |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg_covid <- bind_rows(df_agg_proned_covid, df_agg_unproned_covid) |>
#Only have the propensity at the period and outcome level so need to add period back here
mutate(study_period=fcase(
period_sarscov2=='COVID_SarsCov2Neg-Unk', 'COVID',
period_sarscov2=='COVID_SarsCov2Pos', 'COVID',
period_sarscov2=='Post-COVID_SarsCov2Neg-Unk', 'Post-COVID',
period_sarscov2=='Post-COVID_SarsCov2Pos', 'Post-COVID',
period_sarscov2=='Pre-COVID', 'Pre-COVID'
))
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
#Merge Back with DF Agg COVID
df_agg_covid <- left_join(df_agg_covid, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, period_sarscov2) |>
arrange(hospital_id, period_sarscov2, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, prone_12hour_outcome, study_period)`set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PreMCOVID'),
#COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2Pos'),
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2NegMUnk'),
#Post-COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2Pos'),
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2NegMUnk'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_covid <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type +
(0 + period_sarscov2 | hospital_id)),
data = df_agg_covid,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start samplingplot(agg_brm_period_covid)
summary(agg_brm_period_covid)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type + (0 + period_sarscov2 | hospital_id)
## Data: df_agg_covid (Number of observations: 318)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 37)
## Estimate
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.94
## sd(period_sarscov2COVID_SarsCov2Pos) 0.91
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.97
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.60
## sd(period_sarscov2PreMCOVID) 0.96
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.86
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.79
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.69
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.55
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.52
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.61
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.47
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.56
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.43
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.39
## Est.Error
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.15
## sd(period_sarscov2COVID_SarsCov2Pos) 0.12
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.16
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.29
## sd(period_sarscov2PreMCOVID) 0.24
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.09
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.11
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.13
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.28
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.27
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.27
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.21
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.19
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.22
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.31
## l-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.68
## sd(period_sarscov2COVID_SarsCov2Pos) 0.70
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.69
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.08
## sd(period_sarscov2PreMCOVID) 0.57
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.64
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.52
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.40
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.13
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) -0.16
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.09
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.01
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.13
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) -0.03
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) -0.29
## u-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.27
## sd(period_sarscov2COVID_SarsCov2Pos) 1.18
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.32
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.22
## sd(period_sarscov2PreMCOVID) 1.51
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.97
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.95
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.90
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.93
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.90
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.94
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.85
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.88
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.82
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.87
## Rhat
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2COVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## Bulk_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1989
## sd(period_sarscov2COVID_SarsCov2Pos) 1643
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2405
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1551
## sd(period_sarscov2PreMCOVID) 2325
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1675
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1705
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2180
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 3063
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 3494
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2935
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2101
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 2810
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2733
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1380
## Tail_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 2723
## sd(period_sarscov2COVID_SarsCov2Pos) 2573
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 3002
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1132
## sd(period_sarscov2PreMCOVID) 2778
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 2525
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2774
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2868
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2418
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 2430
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2294
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2790
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 3124
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2749
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1911
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## period_sarscov2COVID_SarsCov2NegMUnk -0.98 0.23 -1.43 -0.53
## period_sarscov2COVID_SarsCov2Pos 0.47 0.22 0.01 0.89
## period_sarscov2PostMCOVID_SarsCov2NegMUnk -1.24 0.24 -1.71 -0.78
## period_sarscov2PostMCOVID_SarsCov2Pos -0.27 0.27 -0.82 0.26
## period_sarscov2PreMCOVID -2.28 0.29 -2.89 -1.74
## Hospital_TypeCommunity 0.83 0.26 0.32 1.33
## Rhat Bulk_ESS Tail_ESS
## period_sarscov2COVID_SarsCov2NegMUnk 1.00 843 1671
## period_sarscov2COVID_SarsCov2Pos 1.00 1012 1929
## period_sarscov2PostMCOVID_SarsCov2NegMUnk 1.00 1032 1737
## period_sarscov2PostMCOVID_SarsCov2Pos 1.00 1878 2574
## period_sarscov2PreMCOVID 1.00 1473 2179
## Hospital_TypeCommunity 1.00 1212 1775
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).#Compare Differences in COVID vs Not COVID in COVID and POST COVID Period
post_covid <- as_draws_df(agg_brm_period_covid)
# Compute OR for the contrast
delta_covid_period <- post_covid$b_period_sarscov2COVID_SarsCov2Pos - post_covid$b_period_sarscov2COVID_SarsCov2NegMUnk
or <- exp(delta_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period \n', posterior_summary)## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period
## 4.240125 3.195933 5.605969# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2Pos - post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period
## 2.657148 1.55299 4.502385#COVID Negative in COVID Period vs Pre-COVID
# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk - post_covid$b_period_sarscov2PreMCOVID
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID
## 2.794761 1.608803 5.018258#ICU bed capacity Not Included as an additional covariable (in addition to hospital type becaues of co-linearity)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Size
prior(normal(0,1), class = "b", coef = 'icu_bed_catlarge'),
prior(normal(0,1), class = "b", coef = 'icu_bed_catmedium'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_size <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + icu_bed_cat +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_size)
summary(agg_brm_period_size)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + icu_bed_cat + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.94 0.13 0.72
## sd(study_periodPostMCOVID) 1.04 0.18 0.73
## sd(study_periodPreMCOVID) 1.02 0.25 0.61
## cor(study_periodCOVID,study_periodPostMCOVID) 0.80 0.10 0.57
## cor(study_periodCOVID,study_periodPreMCOVID) 0.63 0.19 0.20
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.55 0.21 0.09
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.24 1.00 1133
## sd(study_periodPostMCOVID) 1.44 1.00 1494
## sd(study_periodPreMCOVID) 1.59 1.00 1536
## cor(study_periodCOVID,study_periodPostMCOVID) 0.94 1.00 1348
## cor(study_periodCOVID,study_periodPreMCOVID) 0.92 1.00 1848
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.90 1.00 1450
## Tail_ESS
## sd(study_periodCOVID) 2086
## sd(study_periodPostMCOVID) 2381
## sd(study_periodPreMCOVID) 2456
## cor(study_periodCOVID,study_periodPostMCOVID) 2211
## cor(study_periodCOVID,study_periodPreMCOVID) 2499
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1850
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.69 0.25 0.19 1.18 1.00 1106
## study_periodPostMCOVID -0.33 0.28 -0.89 0.20 1.00 1268
## study_periodPreMCOVID -1.41 0.33 -2.07 -0.80 1.00 1409
## icu_bed_catmedium -0.34 0.36 -1.04 0.38 1.01 1022
## icu_bed_catlarge -0.55 0.33 -1.20 0.11 1.00 956
## Tail_ESS
## study_periodCOVID 1840
## study_periodPostMCOVID 1842
## study_periodPreMCOVID 2475
## icu_bed_catmedium 1332
## icu_bed_catlarge 1519
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.496 seconds (Warm-up)
## Chain 2: 2.271 seconds (Sampling)
## Chain 2: 4.767 seconds (Total)
## Chain 2:
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.663 seconds (Warm-up)
## Chain 4: 2.127 seconds (Sampling)
## Chain 4: 4.79 seconds (Total)
## Chain 4:
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.442 seconds (Warm-up)
## Chain 3: 2.425 seconds (Sampling)
## Chain 3: 4.867 seconds (Total)
## Chain 3:
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.887 seconds (Warm-up)
## Chain 1: 2.973 seconds (Sampling)
## Chain 1: 5.86 seconds (Total)
## Chain 1:#Now Generate Odds Ratios for Hospital Size Contrasts
full_posterior <- as_draws_df(agg_brm_period_size) |>
mutate(medium_v_small = exp(b_icu_bed_catmedium ),
large_v_small = exp(b_icu_bed_catlarge))
summary_posterior_size <- data.frame(
variable = c(
'Medium vs Small ICU Capacity',
'Large vs Small ICU Capacity'
),
median_or = c(median(full_posterior$medium_v_small),
median(full_posterior$large_v_small)),
conf.low = c(quantile(full_posterior$medium_v_small, p=0.025),
quantile(full_posterior$large_v_small, p=0.025)),
conf.high = c(quantile(full_posterior$medium_v_small, p=0.975),
quantile(full_posterior$large_v_small, p=0.975)))
print(summary_posterior_size)
write.csv(summary_posterior_size, paste0(project_location, '/summary_output/summary_posterior_size_analysis.csv'))#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_strat <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_strat)
summary(agg_brm_period_strat)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.97 0.13 0.76
## sd(study_periodPostMCOVID) 1.08 0.17 0.79
## sd(study_periodPreMCOVID) 1.01 0.24 0.61
## cor(study_periodCOVID,study_periodPostMCOVID) 0.82 0.09 0.60
## cor(study_periodCOVID,study_periodPreMCOVID) 0.64 0.19 0.21
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.57 0.22 0.08
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.26 1.00 991
## sd(study_periodPostMCOVID) 1.47 1.00 1619
## sd(study_periodPreMCOVID) 1.55 1.00 1949
## cor(study_periodCOVID,study_periodPostMCOVID) 0.94 1.00 1731
## cor(study_periodCOVID,study_periodPreMCOVID) 0.94 1.00 1826
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.91 1.00 1415
## Tail_ESS
## sd(study_periodCOVID) 1675
## sd(study_periodPostMCOVID) 2482
## sd(study_periodPreMCOVID) 3092
## cor(study_periodCOVID,study_periodPostMCOVID) 2655
## cor(study_periodCOVID,study_periodPreMCOVID) 2238
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1708
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.41 0.17 0.07 0.74 1.01 565
## study_periodPostMCOVID -0.62 0.21 -1.04 -0.22 1.00 803
## study_periodPreMCOVID -1.70 0.25 -2.22 -1.23 1.00 1492
## Tail_ESS
## study_periodCOVID 953
## study_periodPostMCOVID 1715
## study_periodPreMCOVID 2623
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## tion: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 1.486 seconds (Warm-up)
## Chain 4: 1.163 seconds (Sampling)
## Chain 4: 2.649 seconds (Total)
## Chain 4:
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.506 seconds (Warm-up)
## Chain 3: 1.166 seconds (Sampling)
## Chain 3: 2.672 seconds (Total)
## Chain 3:
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.787 seconds (Warm-up)
## Chain 1: 1.163 seconds (Sampling)
## Chain 1: 2.95 seconds (Total)
## Chain 1:
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.838 seconds (Warm-up)
## Chain 2: 1.173 seconds (Sampling)
## Chain 2: 3.011 seconds (Total)
## Chain 2:#Posterior Draws
sensitivity_posterior <- as_draws_df(agg_brm_period_strat) |>
mutate(covid_v_pre_covid=exp(b_study_periodCOVID-b_study_periodPreMCOVID),
covid_v_post_covid=exp(b_study_periodCOVID-b_study_periodPostMCOVID),
pre_covid_v_covid=exp(b_study_periodPreMCOVID-b_study_periodCOVID),
post_covid_v_covid=exp(b_study_periodPostMCOVID-b_study_periodCOVID),
post_covid_v_precovid=exp(b_study_periodPostMCOVID-b_study_periodPreMCOVID))
sensitivity_table <- data.frame(
variable = c('COVID_v_Pre-COVID',
'COVID_v_Post-COVID',
'Pre-COVID_v_COVID',
'Post-COVID_v_COVID',
'Post-COVID_v_Pre-COVID'),
or = c(median(sensitivity_posterior$covid_v_pre_covid),
median(sensitivity_posterior$covid_v_post_covid),
median(sensitivity_posterior$pre_covid_v_covid),
median(sensitivity_posterior$post_covid_v_covid),
median(sensitivity_posterior$post_covid_v_precovid)),
conf.low=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.025),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.025),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.025),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.025),
quantile(sensitivity_posterior$post_covid_v_precovid, p=0.025)),
conf.high=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.975),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.975),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.975),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.975),
quantile(sensitivity_posterior$post_covid_v_precovid, p=0.975)))
sensitivity_table
write.csv(sensitivity_table, paste0(project_location, '/summary_output/sensitivity_table.csv'))